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Abstract 

Background: Computed Tomography (CT) is a technology that obtains the 
tomogram of the observed objects. In real-world applications, especially the 
biomedical applications, lower radiation dose have been constantly pursued. To 
shorten scanning time and reduce radiation dose, one can decrease X-ray exposure 
time at each projection view or decrease the number of projections. Until quite 
recently, the traditional filtered back projection (FBP) method has been commonly 
exploited in CT image reconstruction. Applying the FBP method requires using a 
large amount of projection data. Especially when the exposure speed is limited by 
the mechanical characteristic of the imaging facilities, using FBP method may 
prolong scanning time and cumulate with a high dose of radiation consequently 
damaging the biological specimens. 

Methods: In this paper, we present a compressed sensing-based (CS-based) iterative 
algorithm for CT reconstruction. The algorithm minimizes the / 7 -norm of the sparse 
image as the constraint factor for the iteration procedure. With this method, we can 
reconstruct images from substantially reduced projection data and reduce the 
impact of artifacts introduced into the CT reconstructed image by insufficient 
projection information. 

Results: To validate and evaluate the performance of this CS-base iterative algorithm, 
we carried out quantitative evaluation studies in imaging of both software Shepp- 
Logan phantom and real polystyrene sample. The former is completely absorption 
based and the later is imaged in phase contrast. The results show that the CS-based 
iterative algorithm can yield images with quality comparable to that obtained with 
existing FBP and traditional algebraic reconstruction technique (ART) algorithms. 

Discussion: Compared with the common reconstruction from 180 projection 
images, this algorithm completes CT reconstruction from only 60 projection images, 
cuts the scan time, and maintains the acceptable quality of the reconstructed 
images. 
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Background 

Computed Tomography (CT), which obtains a series of projection data of objects con- 
cerned from several view angles, can get the tomograms of the objects through the 
technology of image reconstruction. From a purely mathematical standpoint, the solu- 
tion to the problem of how to reconstruct a function from its projections dated back 
to the paper by Radon in 1917. The current excitement in tomographic imaging origi- 
nated with Hounsfield's invention of the X-ray computed tomographic scanner for 
which he received a Nobel Prize in 1972 [1]. The algorithms of image reconstruction 
from projections can be divided into two classes: the analytical method and the alge- 
braic method [2]. The advantages of the analytical method, such as filtered back pro- 
jection (FBP) method, are relatively high computational speed and short computational 
time. When the projection data are densely sampled, images can be reconstructed 
accurately with analytic methods [3]. Thus, this method is widely used in the commer- 
cial CT systems. However, if data containing a reduced number of projections sparsely 
sampled over an angular range are considered, the analytic algorithms will yield recon- 
structed images with severe aliasing artifacts such as sharp streaks [4]. The iterative 
algorithms, on the contrary, can reconstruct images from relatively less projection data. 
But, it will take much longer time with iterative algorithms versus analytic algorithms. 

In real-world applications, especially the biomedical applications, higher temporal 
resolution and lower radiation dose have been constantly pursued. One can reduce 
scanning time and radiation dose by decreasing X-ray exposure time at each projection 
view. However, the reduction of exposure time would further lower the signal-to- noise 
ratio (SNR) of the projection images and consequently lower the reconstructed images' 
quality [5]. The other approach to decrease scanning time and radiation dose is to 
reduce the number of projections. 

Compressed sensing theory, also known as compressive sampling or CS, suggested by 
Candes, Romberg, Tao and Donoho in 2006 [6,7], states that one can reconstruct 
images accurately from a number of samples that are far smaller than the desired reso- 
lution of the images [8]. Inspired by the theory's success in signal recovery, we have 
anticipated that a CS-based algorithm may be used to reconstruct images from sub- 
stantially reduced projection data. The algorithm minimizes the ^-norm of the sparse 
image as the constraint factor for the iteration procedure. This work focuses on recon- 
structing images from significantly reduced projection data, shortening scanning time, 
minimizing radiation dose without reducing image quality, and employing this algo- 
rithm in phase contrast imaging experiments. 

The paper is organized as follows. In section 2, the experimental set-up will be intro- 
duced for our polystyrene sample imaging. In section 3, the materials and methods for 
data scanning and image reconstruction will be described, in section 4, numerical 
results under different conditions are reported and the reconstructed and reference 
images at the visualization-based and quantitative-metric-based evaluation levels are 
compared. Finally, the implication of the results will be further discussed in section 5. 

The experimental set-up for phase contrast imaging 

Phase contrast X-ray imaging [9-13] enables the observation of light samples, such as 
biological soft tissue, without a contrast agent, because the phase shift cross sections of 
light elements are much larger than their absorption cross sections [14]. However, in 
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the research of phase contrast imaging, the FBP method has been commonly exploited 
in CT reconstruction [15-17]. Applying the FBP method requires a large amount of 
projection data, which can prolong scanning time and cumulate with a high dose of 
radiation potentially damaging the biological specimens. 

Our experiment was performed at the 4W1A beamline of Beijing Synchrotron Radia- 
tion Facilities (BSRF). The synchrotron X-ray beam was 20 mm in width by 10 mm in 
height. The X-ray beam energy was set at 15keV in this experiment. The distance 
between the synchrotron radiation source and the sample was approximately 43 m. 
The X-ray CCD was the Photonic Science X-ray FDI-18 mm camera system with 1300 
x 1030 pixels, and 10.9 x 10.9 urn 2 per pixel. The schematic set-up of this analyzer 
based imaging (ABI) system is shown in Figure 1. It composes of a monochromator 
crystal, an analyzer crystal, one sample rotation stage, and one image detector. The 
monochromator crystal is used to produce highly parallel and monochromatic X-ray 
beams. When the highly parallel and monochromatic X-ray beams travel through the 
object, their directions change due to refraction and scattering. According to the Brag 
diffraction theory, the analyzer crystal only reflects photons coming from a particular 
angle. Thus, if the analyzer crystal is rotated about an axis perpendicular to the meri- 
dian plane, the diffracted intensity will trace out a 'rocking curve' [18]. We obtain the 
projection data of our polystyrene phantom when the analyzer crystal was set to the 
half intensity points on the high-angle sides of the rocking curve. 

Materials and methods 

Both software phantom and real sample were used to test our algorithm. The software 
phantom was a discrete 256 x 256 Shepp-Logan phantom (see Figure 2(a)). We gener- 
ated Sheep-Logan phantom using Matlab command sentence "phantom(256)". We sup- 
pose that it is the desired CT image and each pixel value presents an attenuation 
coefficient. To generate the projection data, we simulated the procedure of an X-ray 
scan, computed line integrals across the image, and eventually, obtained the projection 
data as shown in Figure 3. More popularly, this projection data is named 'sinogram' 



sample 





image detector 



synchrotron X-ray 



rotation stage 



monochromator crystal 



Figure 1 Schematic set-up of an analyzer based X-ray imaging system. 
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(a) (b) 

Figure 2 Shepp-Logan phantom, (a) the original image, (b) the gradient counterpart of (a). 



[19]. The horizontal and vertical axes represent the detector-bin and view- angle coor- 
dinates. In our sample, the number of the detector-bin was 256, and the number of 
the view-angle had two selections which were 60 and 30. Since minimal levels of noise 
were introduced to assess the stability of the algorithm, such studies were supposed to 
give an upper bound to the performance of the CT reconstruction algorithms. The real 
sample was a polystyrene hexahedron. The length, width, and height of this hexahe- 
dron were about 4, 4, and 2 mm respectively. On one 4 mm by 4 mm surface, several 
lines of holes were punched by a laser gun. We placed the other 4 mm by 4 mm sur- 
face on the sample stage. By scanning the sample over 180° with the angular step of 1°, 
we obtained one hundred and eighty projection images of 1300 x 1030 pixels. Several 
pieces of projection images in different angles of rotation are shown in Figure 4. Before 
reconstruction, we moved the background of the projection images. And since the 
sample was not placed at the absolute center of the sample rotation stage, the projec- 
tion images should be cut to a suitable size to compose 'sinogram'. 



I 

Figure 3 Projection data of Shepp-Logan phantom in the 2D data space with the horizontal and 
vertical axes representing the detector-bin and view-angle coordinates, respectively. 
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Figure 4 The Radiographs collected from the Diffraction enhanced X-ray imaging system in 
different rotate angles. 



Now, we describe our CS -based iterative algorithm for image reconstruction in paral- 
lel-beam CT. A successful application of CS requires that the desired image should 
have a sparse representation in a known transform domain [7]. Consider an image f, 
which can be viewed as an N x 1 column vector in R* 1 , whose individual elements fp j 
- 1, 2, ... N are N pixel values of the image. Expand vector fin an orthonormal basis 
¥ as follows: 



(1) 



where ¥ is the N by N matrix [y/±..., y/jq] with the N x 1 vectors {t^-}^ as columns 
and x is also an N x 1 column vector. If all but a few of entries in vector x are zero or 
almost zero, we will say that /is sparse in the ^domain and x is its sparse representa- 
tion. For example, the Shepp-Logan phantom in Figure 2(a) and its gradient counter- 
part in 2(b), we denote the intensity of pixel of a 2D image as f h> w , where h = 1, 2... H; 
w = 1, 2... W; H and W are the height and width of the 2D image respectively and W x 
H = N. If the pixel values are labeled by f h> w , the gradient modulus is as follows. 



^fh,w\ = y (fh+l,w -fh.w) + (fh,w+l -fh,w) 



(2) 



We refer to this quantity as the gradient image [20]. The number of non-zero pixels in 
this 256 x 256 image (Figure 2(a)) is 27521, while the number of non-zero pixels in its 
gradient image (Figure 2(b)) is only 2182, which is much less than the pixel number of 
the image. That means clearly, /and x are equivalent representations of the image, with 
fin the space domain and x in the ^domain. In realistic CT imaging, suppose the 
sampled parallel-beam projection-data of image /are modeled by a discrete linear system 



(3) 



where vector g has length M with individual measurements referred to as g b i = 1, 2, 
... M and 0 is the M by N system matrix [21] that yields the discrete set of projection 
measurements for parallel-beam scanning from the object. Then substituting Wx for/, 
g can be written as 



(4) 
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where ®' = Oij/ is an M by N matrix. For a sparse image, since M < N in Eq. (4) 
there are infinitely many 3c that satisfy g = O/x. Therefore, the image reconstruction is 
aimed at finding the vector x in transform domain by solving the linear program [6-8] 

x = argmin \\x\\ ll subject to &x = g. (5) 

X 

Define the lj -norm of the vector x as = y^ N ^ |x;|[8]. In this paper, specifically, 

x represents the gradient image vector of the desired image. 

The constraints 0'3c = g can be satisfied under the circumstance that measurements 
contain no noise. But it is unattainable, so we consider the constraints that 

\&x — g\ < s, (6) 
where s is a small error factor. 

To minimize the lj -norm of the gradient image [20,22,23], a basic gradient descent 
method was employed. Usually, the gradient descent method is to reduce the objective 
function y 1 = ||g — 0/|| 2 by iteratively moving the image along the gradient [24] 

jrnext _ jrcurrent ^ ^current 

where a is constant to control the descent speed and a is related to the gradient of 
the lj -norm of the gradient image which can also be thought of an image. Each pixel 
value of it is expressed as the following partial derivative [20] 

_ 9 l V /kK'|i 1 _ 2fh,w - fh+l t w ~ fh,w+l 



d f ]hW f e + (fh+hw -fh,u>Y + (fh,w+i -fh,u>Y 



fh,w fh—l,w 

(8) 



£ + (fh,w —fh-hw) + (fh-hw+l —fh-hw) 
fh,w ~ fh,w—\ 

^ + (fh+hw-l — fh,w-l) + (fh,w — fh,w-l) 

To avoid the condition that the denominator vanishes, a small positive number s is 
added in each radical sign. 

In the discrete setting, the parallel-beam projection-data vector g can be written as 
weighted sums over the pixels traversed by the X-ray as 

N 

ft = Z^J where x = 1,2, • • • ,M. (9) 

i=i 

The weight component (\> iy j of the system matrix O is computed by the intersection 
length of the ith ray through the ;th pixel. Using a sketch we can understand it clearly. 
In Figure 5, the image is composed of four pixels f 0 , f lf f 2 , and fy g 1 is a measurement; 
and the X-ray l x passed pixels fo,f\, and^; the lengths passed these three pixels are 0 1} 
o, 0i, i and 0i, 3 respectively. The computation of the weight 0 is complex. The 
immediate computation of each 0 will prolong the reconstruction time. Especially with 
iteration times' increasing, the computation of weights will repeat again and again. A 
solution is to save the weights in a file in advance and read the weights from this file 
during the iteration. Since the X-rays are parallel-beam, we can take advantage of the 
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/b 





Figure 5 Sensing matrix model. f 0/ f h f 2 , and f 3 : four pixels. 0 7O/ 0 77 , and (p 13 : three weight coefficients. 
/; is one X-ray, and g? is one projection value. 



symmetrical characteristic of X-rays. In Figure 6, the X-rays are denoted by a, b, c, and 
d. Actually, the measurements obtained by the integral along the X-rays a h b h Cx and 
dx are detected by the same detector in different rotate angles which are a, (90- a), 
(90+ a) and (180- a) degrees respectively. The a h a 2l ^ a m is a set of parallel beams. 
These eight X-rays [a h a m , b h b m) c h c m , dj and d m ) in Figure 6 have such symmetri- 
cal characteristics that line y - -x is the symmetrical axis of ax and bj; the x-axis is the 
symmetrical axis of bj and c 1} also a ± and df, and the parallel beams a ± and a m , b 1 
and b m , Cx and c m , dj and d m , are symmetrical to the origin. So, if we have a weight 
value in X-ray a 2 (the wide segment in line aj), we can gain the other seven weight 
values (the wide segment in the other seven lines) using the symmetrical 
characteristics. 

In the following, we give the steps of the reconstruction algorithm in the form of a 
pseudo-code and abbreviated notation. 

(1) initialization of the image f: 

/(°) = 0; 

(2) iterative process: 



■(fe) 



N 

J2 <t>i,n • fn 
n=l 



(fe-1) 



N 
i=l 



where the relaxation parameter X [25] is a positive real number to adjust the itera- 
tive, and k is from 1 to M. When k = M , a complete iteration period was finished. 
The next iteration will enforce the estimated image to the constraint — g\ < s by 
the gradient descent iteration. 
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Figure 6 Computation of sensing matrix weights model. Eight arrow lines represent X-rays. The rays at 
different angles are labeled by different letters a, b, c, and d. The subscripts represent the detector array's 
position, a, (90- a), (90+ a) and (180- a): angles from X-axis to X-rays. 

V ) 



(3) initialization of the gradient descent image: 

(4) gradient descent iteration: 

a . a,, and 

A = |/(°) -/ (0) | • Vy,x . 

| Vy,x | 

In this iteration, the end time we selected is 5. 

(5) Initialize the next iterative step: 

y(0) = j{end) i 

then we repeat step (2) - (5) until the difference between the current f M) and the 
previous/^ is smaller than the threshold we set or the iteration number is more than 
1000. 

About the control parameters, we selected A = 1.0, s = 0.0001, and a - 0.5 respec- 
tively. The threshold value to stop iteration was set as 0.001. These presetting para- 
meters and coefficients only appear to alter the convergence rate. 
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Results 

To demonstrate this CS-based iterative algorithm for image reconstruction from 
under-sampled projection data, we performed two sets of studies: the first set of stu- 
dies were designed in such a way as to acquire some theoretical understanding of how 
the CS-based iterative algorithm performs on image reconstruction from reduced pro- 
jection data with the parallel-beam configuration under ideal conditions, and the sec- 
ond set of numerical examples aimed to see how the CS-based iterative algorithm 
could be applied to phase contrast CT image reconstruction. 

Recall Eq. 3, in the situation that the measurement data g contain no noise and the 
full scan views data are used, one might expect to reconstruct images accurately. How- 
ever, in the present studies, the projection data were under-sampled in view angle. We 
performed the FBP, traditional algebraic reconstruction technique (ART) and CS-based 
iterative reconstruction algorithms under the condition that the numbers of views were 
60 and 30. The image-quality evaluation for each specimen was performed at two dif- 
ferent levels, including 1) visualization-based evaluation, and 2) quantitative-metric- 
based evaluation. Some of the evaluation concerns make a comparison between the 
reconstructed and original images. 

Visual inspection of reconstructions in Figure 7 suggests that under the conditions of 
few-view (60- and 30-view number) projection data, the CS-based iterative algorithm 
can effectively suppress streak artifacts and noise observed in images obtained with the 
FBP and traditional ART algorithms, thus yielding images with a higher visual similar- 
ity to the Shepp-Logan phantom image (see Figure 2(a)) than those obtained with 
other algorithms. 

In addition to visualization-based evaluation, the following three metrics were 
employed to quantitatively assess the similarity between reconstructed images and the 
original phantom image: 1) the root mean squared error (RMSE), 2) the universal qual- 
ity index (UQI) [26], and the correlation coefficient (CC), which are defined as 



RMSE : 



N 2 

J2{fn-foi) no) 



2Cov\f ri f 0 } 2f r -f 0 

UQI = V — 1 J L , (11) 

D(f r ) + D(f 0 )f2 +f 2' 
2C0V \fr,fo\ 

CC = 1 , J — , (12) 

where vector f r and f 0 denote the reconstructed and original images of N pixels, and 

I N i N 

/o = — J2foh fr = — J2fri' 

iN/ i=l iN/ i=l 

DCfo) = ^ | (/« -/o) 2 , D(f r ) = ^ | (/„ -l) 2 , 

Cov{f r ,fo} = j ±-J2(f«-~fr) (/«-/»)• 

1=1 
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Figure 7 Shepp-Logan phantom images reconstructed from 60 and 30 view numbers using CS- 
based iterative algorithm (column 1), ART (column 2), and FBP (column 3). 



The RMSE is widely used for measuring reconstruction accuracy, whereas the UQI 
and CC can be used for evaluating the pixel-to-pixel similarity between reconstructed 
and original images. When assessing the image's quality, we demand the RMSE index 
to be as small as possible, while expecting the UQI and CC to have the contrary 
results. 

In Figure 7, the images in the left column are the reconstructed images using CS- 
based iterative algorithm, the middle column using ART, and the right column using 
FBP. The images in row 1 are reconstructed from 60-view data, and in row 2 the 
images are reconstructed from 30-view data. 

From the digital Shepp-Logan phantom and reconstructed images, we computed 
their RMSEs, UQIs, and CCs and summarized them in Table 1, 2, and 3 respectively. 
Results of these three metrics suggest that the CS-based iterative algorithm yields 
images more similar to the original image than the FBP and traditional ART 
algorithms. 

In addition to the simulated experiments with Shepp-Logan phantom, phase contrast 
X-ray imaging of a real polystyrene phantom was also performed. We collected 180 
radiographs at 180 views. From this full data set, we applied the conventional FBP 
algorithm to yield the reference CT image in the 44th slice (see Figure 8). Then we 
extracted the 60- and 30-view subsets of data evenly distributed over n-view to simu- 
late data collected at a reduced number of projection views. The reconstructed images 
are shown in Figure 9, and the three quantitative metrics RMSEs, UQIs, and CCs are 



Table 1 Quantitative assessment RMSE of Shepp-Logan phantom images 



RMSE 


CS-based 


ART 


FBP 




Iteration 






60-view number 


16.82 


22.61 


37.92 


30-view number 


17.26 


25.07 


40.76 
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Table 2 Quantitative assessment UQI of Shepp-Logan phantom images 

UQI CS-based ART FBP 

Iteration 

60-view number 0.942 0.894 0.570 

30-view number 0.938 0.822 0.490 



Table 3 Quantitative assessment CC of Shepp-Logan phantom images 



cc 


CS-based 


ART 


FBP 




Iteration 






60-view number 


0.947 


0.900 


0.719 


30-view number 


0.945 


0.891 


0.646 



computed and summarized in Table 4, 5, and 6 respectively. The results suggest that 
the CS-based iterative algorithm can yield the most similar images to the reference 
image. Different from the results of the reconstructed Shepp-Logan images, images 
reconstructed from the traditional ART algorithm seem to have the least similarity to 
the reference image. The reason might be because the reference image is reconstructed 
from the FBP algorithm as well. 

Discussion and Conclusions 

In this article, a CS-base iterative algorithm reconstructing images from substantially 
reduced projection data was presented. Both the digital Shepp-Logan phantom and the 
real polystyrene sample experiment results show that the CS-based iterative algorithm 
can yield images with quality comparable to that obtained with existing FBP and tradi- 
tional ART algorithms. However, when the number of gradient descent iterations is 
increased, the smoothing artifact in the reconstructed images will be more obvious. To 
improve this situation, reducing the number of gradient descent iterations or the step 




Figure 8 The reference Polystyrene image reconstructed from 180-view projection data using FBP 
algorithm. 
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Table 4 Quantitative assessment RMSE of Polystyrene phantom images 

RMSE CS-based ART FBP 

Iteration 

60-view number 7.47 17.90 12.16 

30-view number 9.52 19.44 17.37 



Table 5 Quantitative assessment UQI of Polystyrene phantom images 



UQI 


CS-based 


ART 


FBP 




Iteration 






60-view number 


0.897 


0.648 


0.797 


30-view number 


0.817 


0.584 


0.613 



size is an alternative, but if the gradient descent is too small, the algorithm will reduce 
to the traditional ART algorithm. 

Because the system matrix for sampling is saved in a matrix file, we must create a 
correlating matrix file for different images in different dimensions at first. The size of 
such a file, however, relates to the projection data's number and images' dimensions. 
For this reason, if either of these two elements is rather big, the matrix file will be 
quite large. Therefore, either compressing the matrix file further or diminishing the 



Table 6 Quantitative assessment CC of Polystyrene phantom images 



cc 


CS-based 


ART 


FBP 




Iteration 






60-view number 


0.900 


0.668 


0.803 


30-view number 


0.831 


0.601 


0.635 
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region of interest (ROI) in the reconstructed images may be one direction of our 
further work. 

In conclusion, the paper aims to reduce the impact of artifacts introduced into the 
CT reconstructed image by insufficient projection information. The feasibility of this 
method which enforces the gradient descent for constraints in traditional iterative algo- 
rithms has been demonstrated by both the simulated phantom and the real polystyrene 
sample experiments. The results show that the CS-based iterative algorithm can yield 
images with quality comparable to that obtained with existing FBP and traditional 
ART algorithms. Further research will be performed to improve algorithm efficiency. 
Moreover, applying this algorithm to a less "sparse" sample such as the real biological 
soft tissues of small animals and studying how effective the method would be is our 
future concern. 
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